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ABSTRACT 

Hydrogen recombination lines are one of the major diagnostics of H II region phys- 
ical properties and kinematics. In the near future, the Expanded Very Large Array 
(EVLA) and the Atacama Large Millimeter Array (ALMA) will allow observers to 
study recombination lines in the radio and sub-mm regime in unprecedented detail. 
In this paper, we study the properties of recombination lines, in particular at ALMA 
wavelengths. We find that such lines will lie in almost every wideband ALMA setup 
and that the line emission will be equally detectable in all bands. Furthermore, we 
present our implementation of hydrogen recombination lines in the adaptive-mesh 
radiative transfer code RADMC-3D. We particularly emphasize the importance of 
non-LTE (local thermodynamical equilibrium) modeling since non-LTE effects can 
drastically affect the line shapes and produce asymmetric line profiles from radially 
symmetric H n regions. We demonstrate how these non-LTE effects can be used as a 
probe of systematic motions (infall & outflow) in the gas. We use RADMC-3D to pro- 
duce synthetic observations of model H II regions and study the necessary conditions 
for observing such asymmetric line profiles with ALMA and EVLA. 



1 INTRODUCTION 



Massive stars influence their cosmic environment through 
powerful winds, radiation feedback and supernova explo- 
sions. Ionizing radiation from massive stars produces pro- 
nounced H II regions around them. Observations of H II 
regions, while still in the ultracompact phase (with diame- 
ters less than 0.1 pc), provide important insight into massive 



H II region expansion during massive star formation ( Peters 
|et al.|2010a]c|b|p0Tl) . 

These collapse simulations require adaptive-mesh simu- 
lations with 10 refinement levels and more. To post-process 
such data, dedicated radiative transfer tools must be devel- 
oped. To create synthetic recombination line observations 
of the simulated H II regions, we have implemented recom- 
bination lines in RADMC-3D. RADMC-3D is a radiative 



star formation i 


Habing & Israeli 1979 Churchwell[2002 Zin- 


|necker & Yorke 


2007l. Observations of hydrogen recombina- 



tion lines at radio and sub-mm frequencies (RRLs) are rou- 
tinely used by observers to infer densities, temperatures and 



transfer code that can handle both continuum (e.g. Peters 
et al.pOlQc l and line (e.g. Shetty et al.||20lT i radiation on 



velocity structures inside H II regions ( Gordon & Sorochenko 
[2002 ). 



While the microphysics of recombination lines was rea- 
sonably well understood relatively early ( |Dupree fc Gold^| 
berg j [l970| l, the H II regions themselves had to be modeled 



arbitrary octree meshes. It has an interface for PARAMESH 
( MacNeice et al.|2000 I , the grid library of the adaptive-mesh 
code FLASH flFryxell et al.|2000| >, which was used for these 
simulations. We have used RADMC-3D previously to model 
free-free and dust continuum emission from the simulated 



either on scales larger than the scale of gravitational col- 
lapse using numerical radiation-hydrodynamics (see reviews 



H ii regions ( jPeters et al.|2010c l. 

To the best of the authors' knowledge, the only other 
codes that can model recombination line observations are 
MOLLIE ( |Keto|1990| >, which we have used in previous work 



by |Yorke||1986| |Mac Low||2008| |Klessen et al.||2011[) or on 

smaller scales with simple spherically JBrown et al. 1978 



to post- process simulation data ( Peters et al 



radiative transfer modeling (Longmore et al 



[Keto|[2002l [2003] ) or non-spherically (Koto 2007]~symmet- 
ric models. Only recently, collapse simulations of massive 
star formation with ionization feedback have left the con- 
striction of two dimensions (Richling & Yorke 1997J and 
facilitated fully three-dimensional dynamical simulations of 



innominate code that |Martm-Pintado et al.| ( 1993| used for 



2010a[ > and for 
2011b, and an 



recombination line models of MWC 349. The advantages of 
RADMC-3D are its direct compatibility with many major 
simulation codes and its modularity. 

Another important feature of RADMC-3D is the pos- 
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sibility to create user-defined setups. Thus, it can be used 
by observers interested in creating simple analytic models of 
the regions they are observing. In the present paper, we de- 
scribe the implementation of hydrogen recombination lines 
in RADMC-3D and demonstrate its capability with syn- 
thetic observations of several H n region models. The devel- 
opment of such a tool is particularly timely because of the 
commissioning of the Expanded Very Large Array (EVLA) 
and the Atacama Large Millimeter Array (ALMA). These 
facilities are set to open new frontiers in terms of sensitiv- 
ity, angular resolution, dynamic range and image fidelity in 
the cm to sub-mm wavelength regimes. Tools that can help 
to interpret such new observational data are clearly highly 
desirable. 

The purpose of the present paper is twofold. First, we 
present our implementation of hydrogen recombination lines 
in RADMC-3D and describe a suite of tests of this imple- 
mentation. The user-defined analytical setups are ideal for 
code validation before applying the method to the much 
more complex simulation data. Second, we use these ana- 
lytical models to simulate synthetic EVLA and ALMA ob- 
servations. We particularly discuss the appearance of line 
asymmetries expected for ALMA observations. 

The paper is organized as follows. Section [2] briefly sum- 
marizes the physics of hydrogen recombination lines. In Sec- 
tion [3] we describe the implementation of recombination 
lines in RADMC-3D. In Section |4) we discuss the general 
properties of RRL transitions at ALMA frequencies. We 
then focus on the observation that the RRL profiles can 
be asymmetric (Section [5| and investigate the conditions 
under which we can expect to observe such asymmetries 
(Section Rjl. We conclude our paper with a summary of our 
results in Section [7] Additional information on line profile 
asymmetries can be found in Appendix [X] We present our 
code tests in Appendix |B| 



2 PHYSICS OF HYDROGEN 
RECOMBINATION LINES 

We start with briefly summarizing the physics of free- 
free continuum radiation and hydrogen recombination lines 



in the radio and sub-millimeter regime (see Gordon & 
Sorochenko (20021 for more details). The hot plasma in 
an H II region gives rise to the emission of thermal 
bremsstrahlung. This free-free radiation causes a continuum 
opacitjQat frequency v of 



a„,c = 0.212 



( — — — \ ( ni ,) 
V 1 cm 3 / V 1 cm -3 / 



T 



(i) 



v 

1 Hz 



with the electron and ion number densities n c and ni, re- 
spectively, and the gas temperature T. Likewise, the plasma 
emits radiation with an emissivity of 



jV,c = B v {T)a v , c , 



(2) 



where B V (T) denotes the intensity of a blackbody of tem- 
perature T at frequency v. 

1 In the following, the subscript C stands for "continuum" , 
whereas the subscript L means "line" . 



During recombination, an excited hydrogen atom emits 
recombination line radiation when the electron falls from at 
state with principal quantum number m to a state with prin- 
cipal quantum number n. The corresponding line absorption 
coefficient is 

) (3) 



47T 



with the Planck constant h, the line profile function <j> v , 
the Einstein coefficients B„, m and Bm,n for absorption and 
stimulated emission, respectively, and the number densities 
of atoms iVfc in state k (for k — m,n). The line emissivity is 

hv , . 

>.l = -7— 4>uJS m A mn (4) 

47T 

with the Einstein coefficient for spontaneous emission 

.3 

(5) 



A - 2hv B 



with the speed of light c. 

The Einstein coefficients satisfy the relation 

giBi^ = gkB k ,i (6) 

with the statistical weights g k = 2k 2 . They are usually ex- 
pressed in terms of the oscillator strengths 



_ m e chy 

471"^ 



(7) 



with the electron mass m e and the electron charge e. The 
oscillator strengths for hydrogen can be approximated for 
large n and small An = m — n following Menzel ( 1968 1 as 



fn,r> 



nM An 1 + 1.5 



.An 



(8) 



with Mau = 0.190775, 0.026332, 0.0081056, 0.0034917, 
0.0018119, 0.0010585 for An = 1,2,3,4,5,6, respectively. 

The situation is more complicated for recombination 
lines in the ALMA bands because lines above 100 GHz show 



non-negligible fine-structure splitting (Towle et al. 19961 



For each principle quantum number k there are k different 
energies which depend on the total angular momentum j — 
1/2, 3/2, . . . ,k — 1/2. These k energies have a degeneracy 



9k,i 



2(2.7 + 1) jsCfc-3/2 



2j + l 



j = k- 1/2. 



(9) 



The fine-structure components typically span a frequency 
range of the order MHz and thus cannot individually be re- 
solved with ALMA due to their thermal broadening in the 
ionized gas, which can be much larger than their separa- 
tion. However, the splitting can change the line profile, and 
the frequencies and Einstein coefficients must be computed 
using relativistic quantum mechanics ( Towle et al.|1 996). 

The line profile function </>„ is in general a convolution 
of a Gaussian profile cf)^ caused by thermal and microtur- 
bulent broadening and a Lorentzian profile 0^ caused by 
electron pressure broadening ( [Brocklehurst fc Seaton|1972[ ). 
The Gaussian profile is given by 



'2-KCF 



exp 



2(7 2 



with 



2c 2 



1 'ZksT , a2 



(10) 



(11) 
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where m p is the proton mass, fee the Boltzmann constant, 
fo the rest frequency of the line and £ the root-mean-square 
of the microturbulent velocity field. The Lorentzian profile 
can be written as 



(12) 



7T (y — vo) 2 + 5 2 



The scale parameter 8 can be approximated following |Brock-| 
lehurst & Leeman ( 1971 1 as 



4.7 



(— V 

V 100/ 



T 



10 4 K 



V 1 cm" 3 ) 



Hz. 



(13) 



The convolution of these two profiles then gives a Voigt pro- 
file ( [Rybicki fc Lightman|1979"t 

1 -H(a,x) (14) 



with a — 8/ y2a, x = (y 
H(a, x) = - 



2tto- 

vo)/\/2o and the Voigt function 
t 2 )dt 



exp(. 



a 2 + (t - x) 2 



(15) 



To determine the occupation numbers, we first act on 
the assumption of local thermodynamic equilibrium (LTE) 
and then consider possible departures from the LTE occu- 
pation numbers. Under the assumption of LTE, the number 
density of hydrogen atoms in state k with energy E k is given 
by the Saha-Boltzmann equation 

k 2 h 3 



Nt 



n c rii 



' T 3 / 2 (27rmfc B ) 3 / 2 6XP \ksT , 
Here, we follow the the sign convention of Gordon fc| 



E k 



(16) 



Sorochenko 



(20021, where Eh > 0, and the factor k derives 
from the statistical weight g k . See Osterbrock (1989) for a 
derivation of this equation. With these occupation number 
densities, we can evaluate the absorption coefficient Q and 
the emissivity 



• LTE D ,rr\ LTE 



(17) 



In general, the occupation numbers will be different from 
the LTE values. This deviation can be measured with a de- 
parture coefficient 



N k = b k Nt- 



(18) 



The coefficient b k represents the fractional departure of Nk 
from A^ TE , with b k = 1 meaning no deviation at all. From 
equation Q it follows that the emissivity under non-LTE 
conditions is 

>,l = bmjv™, (19) 

whereas the absorption coefficient Q is usually written in 
the form 



Q„,l 



b n p n 



with 



Pn,n 



exp 



hp 
knf 



exp 



hv \ 



(20) 



(21) 



In the case of fine-structure splitting, we assume that the g k 
states g k j that correspond to principal quantum number k 



are equally occupiecj^] Thus, the occupation number for the 
fine-structure component with total angular momentum j is 
Nk,j = g k ,jN k /g k . 

The radiative transfer problem for recombination lines 
is commonly solved under the approximation that Thomp- 
son scattering can be neglected. We thus arrive at the ra- 
diative transfer equation 



ds 



(22) 



for the intensity I v at frequency v. With the optical depth 



T v (r) 



{oiu,c{r') + a v ,h(r')) dr' 



o 



and the source function 



S v = 



Equation ( 22 1 can be written as 



Iu{t v ) 



e T " S v (tI) dr^, 



(23) 



(24) 



(25) 



where we have neglected background irradiation. 



3 IMPLEMENTATION 

We have described our implementation of free-free radia- 
tion in RADMC-3D earlier ( jPeters et aL"1|2010c[ ), so that 
we focus here on the implementation of recombination lines. 



The Voigt function ( 15 1 is calculated with the source code 
provided by |Schreier ( jl992] ). The algorithm is an improved 



version of the rational approximation method originally de- 



veloped by Humlfcek ( 1982 1 



The departure coefficients b k for the occupation num- 
ber density N k of quantum number k axe pre-calculated 
in tabulated form with a program published in Gordon & 
Sorochenko ( 2002 I . Since the coefficients vary smoothly with 
logT and logn e , they can be interpolated bi-linearly from 
the tabulated values during the raytracing. The algorithm 
to calculate the b k is based on the matrix condensation tech- 



nique that was originally presented in |Brocklehurst ( 1970 \ 
and Brocklehurst & Salem ( 1977 1 and later extended by 



Walmsley ( 1990 1 towards smaller quantum numbers. In ad- 
dition to the dependence on T and n e , b k also depends on 



the ambient radiation field (see Brocklehurst & Salem 19771 



The radiation field is characterized by the radiation temper- 
ature (10 4 K) and the emission measure (10 6 pccm -6 ). 

We have verified the program by reproducing some 



of the numerical values tabulated in Brocklehurst (19701, 



and 



Brocklehurst & Salem (19771, Salem & Brocklehurst (19791 



Walmsley (19901. As an example, we present in Fig- 



ure[ljthe variation of b k with principal quantum number k at 
T — 10 4 K for different electron number densities n c ranging 

2 We could also assume that the fine-structure levels follow 
a Boltzmann distribution. However, the fine-structure splitting 
changes the energy only by a very small fraction and the dif- 
ference in the energy levels is weak compared to the variation of 
the corresponding Einstein coefficients. Hence, given that all fine- 
structure components overlap due to thermal broadening, we do 
not expect to see noticeable differences in this case. 
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Figure 1. Departure coefficient fej. as function of the principal 
quantum number k at temperature T = 10 4 K for electron number 
densities n e between 10 4 and 10 8 cm - 3 . Compare Figure 1 in 
|Walmsley| | |l990| l. 




50 100 150 200 250 

k 

Figure 2. Departure coefficient bk as function of the principal 
quantum number k at temperature T = 10 4 K for electron num- 
ber densities n c between 10 and 10 4 cm -3 . Compare Figure 2 in 
|Sejnowski fc Hjellmiiig| \1969\ . 



from 10 4 to 10 s cm . The resulting plot is identical to Fig- 



ure 1 in Walmsley (1990). To furthermore estimate the varia- 



tion of the results obtained from different methods, we com- 
pare the bk numbers to values that Sejnowski & Hjcllming 
( 1969 1 computed independently. Figure [2] shows bk as func- 
tion of k at T = 10 4 K for n c between 10 and 10 4 cnT 3 . The 



deviation from Figure 2 in Sejnowski & Hjellming (1969) 



is as small as expected given the different method of com- 
puting the departure coefficients, but according to [Gordon] 
fc Sorochenko| (|2002[) the Walmsley ( 1990 1 values are more 



accurate. 



Gordon & Sorochenko ( 2002 1 argue that the depar- 



ture coefficients calculated by Storey & Hummer ( 1995 \ are 



the most precise ones at the smallest quantum numbers, 
where deviations between the different calculations occur. 



However, the discrepancy between the Walmsley ( 1990 1 and 



Storey & Hummer ( 1995 I coefficients is less than 10% in the 



case they discuss. Given this relatively small difference, we 
will work in this paper with thelWalmsley ( 1990 1 coefficients 



but plan to use more accurate numbers in future work. 

To calculate the fine-structure splitting of hydrogen, 
we have used the code presented in |Gordon fc Sorochenko] 
( 2002 1. We have verified the implementation by reproducing 

19961 for the frequen- 



the values in Table 3 of Towle et al. 



cies and intensities of some recombination lines. 

We refer to Appendix [B] for a presentation of our code 
verification. 



4 GENERAL PROPERTIES OF RRL 

TRANSITIONS AT ALMA FREQUENCIES 

Since the capabilities of ALMA for observing RRL transi- 
tions have not been studied in detail yet, we make some 
general remarks on the properties of recombination line 
observations with ALMA. Once fully operational, ALMA 
will observe at ten bands with wavelengths from 9.7 mm to 
0.3 mm (frequencies of 31 to 950 GHz, see Table [TJ. Note 
that ALMA band 1 overlaps with the upper frequency range 
of the EVLA. The angular resolution in column 4 gives the 



expected range in synthesized beams between the most com- 
pact and most extended array configurations. The line sen- 
sitivity depends on the synthesised beam size (angular res- 
olution), as well as on the receiver sensitivity, weather con- 
ditions (characterized by the precipitable water vapour con- 
tent), the integration time and the spectral resolution. The 
values in column 6 correspond to the sensitivities assuming 
default weather conditions for that frequency (determined 
from the ALMA Observing TooQ, an integration time of 
lmin and a spectral resolution of lkrns" 1 . The two sensi- 
tivies correspond to values for the compact and the most ex- 
tended array configuration. Table entries marked with "*" 
correspond to ALMA bands that are still under develop- 
ment. The data for this table are taken from the ALMA 
Early Science Primei]^] 

To systematically investigate the properties of recom- 
bination lines at ALMA frequencies, we have generated a 
list of all recombination lines with 20 ^ n ^ 200 and 
1 ^ An sC 6. Neglecting stimulated emission, the line 
strength will be determined by the line emissivity (see equa- 
tion Q in Section [2J. Hence, we use the quantity 5*l = 
N m A m , n = A r m 6 m 7V m TE with m = n + An as a measure 
of the line strength. We assume an electron number density 
n c = 5 x 10 4 cm" 3 and a temperature T = 10 4 K inside the 
H II region. For these fiducial values, we calculate St for the 
recombination lines, focusing on the relatively unexplored 
ALMA frequency range. 

It is interesting to note that two effects partially can- 
cel each other in their contribution to the line strength. On 
the one hand, the Einstein coefficient A m>n grows with fre- 
quency because larger frequencies at fixed An correspond to 
smaller quantum numbers m. On the other hand, the states 
with smaller quantum numbers are less occupied under H II 
region conditions than states with higher quantum numbers, 
so the occupation number A^ TE falls off with increasing fre- 
quency. The simultaneous consideration of these two effects 



3 http: / / almascience.eso.org/ call-for-proposals/observing- 
tool/ observing- tool 

4 http: / / almatelescope.ca/ALMA-ESPrimer.pdf 
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Table 1. ALMA full array specifications 



] ) ; \ 1 1 ( | 


wavelength 
(mm) 


frequency 
(GHz) 


angular resolution 
(arcsec) 


continuum sensitivity 
(mjy/beam) 


line sensitivity 
(K) 


1 


6.7-9.5 


31.3-45 


13-0.1 


* 


*/* 


2 


3.3-4.5 


67-90 


6-0.05 


* 


*/* 


3 


2.6-3.6 


84-116 


4.9-0.038 


0.05 


0.07/482 


4 


1.8-2.4 


125-163 


3.3-0.027 


0.06 


0.071/495 


5 


1.4-1.8 


163-211 


* 


* 


*/* 


6 


1.1-1.4 


211-275 


2.0-0.016 


0.10 


0.104/709 


7 


0.8-1.1 


275-373 


1.5-0.012 


0.20 


0.29/1128 


8 


0.6-0.8 


385-500 


1.07-0.009 


0.40 


0.234/1569 


9 


0.4-0.5 


602-720 


0.68-0.006 


0.64 


0.641/4305 


10 


0.3-0.4 


787-950 


0.52-0.005 


1.2 


0.940/* 



Table entries marked with "*" are not yet available. 



order An is (almost) always larger than the highest inten- 
sity line with transition order An + 1 in any single ALMA 
band, but not necessarily if several bands are combined. As 
the figure shows, there are many transitions covering the 
entire ALMA frequency range, so there will likely be a re- 
combination line in any wideband frequency setup. The re- 
combination lines are distributed almost equally with the 
logarithm of frequency. For the convenience of the reader, 
we provide tables containing all lines shown in Figure [4] as 
supplementary online material. 

We note that the true line brightness can only be deter- 
mined by radiative transfer calculations through model H II 
regions that also take stimulated emission and absorption 
into account. We present such calculations in Section |5.1 I 
For one example H II region, |Gordon &; Sorochen ko ( 2002 ) 
calculate the effect of stimulated emission under non-LTE 
conditions and find line depression for ct-transitions in the 
ALMA range, with higher-frequency lines being more de- 
pressed than lower-frequency lines. However, this effect is 
not strong enough to interfere with the general trend shown 
in Figure [4] 




10 2 r 



10 1 f 



1 



10" 



100 200 300 400 500 600 700 
v (GHz) 



Figure 3. The departure coefficient b m , the level population 
Nf^ (scaled by a factor of 10 10 ), Einstein coefficient A m ^ n and 
line strength Sl (scaled by a factor of 10 10 ) for a-transitions in an 
H II region with n E = 5 X 10 4 cm -3 and T = 10 4 K for frequencies 
between 30 and 720 GHz. 



(together with a small decrease of the departure coefficient 
b m with growing frequency) leads to an increase in the line 
strength Sl of more than one order of magnitude from the 
lowest to the highest frequency ALMA bands. As an illus- 
tration, Figure [3] shows the contributions of all three factors 
to Sh across the ALMA frequency coverage. 

As Figure [3] indicates, recombination lines in the higher 
frequency bands will be stronger than for lower frequencies. 
However, as shown in Table [TJ the expected ALMA sensitiv- 
ity will decrease with increasing frequency. Since both the 
line strength and the expected thermal noise grow roughly 
by a factor of ten across the ALMA frequency range, we ex- 
pect the recombination lines of a fixed transition order An 
in all ALMA bands to be equally detectable. 

The rest frequencies of individual transitions in the 
ALMA bands are plotted in Figure [4] We have also marked 
the rest frequencies of some important transitions of the CO 
molecule. For orientation, we show the lower principal quan- 
tum number n for the first and last transitions in a given 
transition order (at the left and right edges of the plot). As 
expected, the lines become brighter as the transition order 
An decreases. The lowest intensity of a line with transition 



5 ASYMMETRIC LINE PROFILES IN MM 
AND SUB-MM RRL SPECTRA 

One interesting feature of hydrogen recombination lines in 
the millimeter and sub-millimeter wavelength regime, as op- 
posed to recombination lines at centimeter wavelengths, is 
the prevalence of non-LTE conditions. Appendix [A] explains 
how this can give rise to asymmetric line profile shapes. 
Essentially, the non-LTE effect is similar to a temperature 
gradient, and the resulting line profiles resemble the well- 
known asymmetric (P-Cyg or inverse P-Cyg) profiles which 
result from optically thick lines in expanding or collaps- 
ing envelopes with a temperature gradient. As the non-LTE 
emission becomes optically thick, the shape of the asymme- 
tries can provide important information about systematic 
motions in the gas. Therefore, it is useful to know which 
physical conditions will cause a given line to i) deviate from 
LTE and ii) become optically thick. 

The deviation of a line from LTE is determined by the 
departure coefficient bk (see Equation (18 1 in Section [2j|. 
Values further from one correspond to larger deviations from 
LTE. Figure[5]shows the departure coefficients bk as function 
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Figure 4. Intensity of hydrogen recombination line transitions in the ALMA bands. The shaded regions mark the (partly overlapping) 
ALMA bands. The dashed lines show the location of some CO transitions. The numbers on the left- and right-hand side indicate the 
lower quantum number n for the first and last transition for a series of fixed An, respectively. 



of frequency for transitions with An = 1 (giving a unique 
relationship between k and 6&) and electron number densi- 
ties between 10 4 and 10 s cm -3 . The shaded regions mark the 
ALMA bands. It becomes clear that departures from LTE at 
typical H II region densities ( Kurtz|2005 I for hypercompact 
(n c > 10 6 cm -3 ), ultracompact (n e > 10 4 cm -3 ) and com- 
pact (n e > 10 3 cm -3 ) regions become larger at frequencies 
above 100 GHz, or wavelengths shorter than 0.3 cm. 

Departure coefficients smaller than unity are a neces- 
sary, but not sufficient condition to see strong non-LTE ef- 
fects. The emission optical depth plays an important role 
too, as optically thin emission will not give rise to line pro- 
file asymmetries. Determining whether strong asymmetries 
are produced requires radiative transfer modelling of the 
RRL emission. In Section [5. 1| we generate simple analytical 
H II regions, solve the non-LTE radiative transfer equations 
for RRLs at frequencies covered by the EVLA and ALMA, 
then simulate observing these regions for typical EVLA and 
ALMA observing setups. 



5.1 Simulated EVLA and ALMA observations 

We aimed to simulate RRL observations of H II regions, 
covering a range of typical EVLA and ALMA observational 
setups. Synthetic spherically symmetric H II regions were 
generated using analytical profiles for n c (equation (|B1[)) 



and the radial velocity v r (equation (B6l), corresponding to 



model 6 from Appendix [BJ which were chosen as represen- 
tative of compact H II regions. We then used the radiative 
transfer modelling to generate synthetic emission maps of 
the H II regions. The radiative transfer models deal with 
continuum and line emission simultaneously so create a com- 
bined line plus continuum emission cube. We separated the 



two components by fitting a low-order polynomial to the 
line-free channels and from there on treated each compo- 
nent individually. In this way synthetic continuum and RRL 
emission cubes were generated for transitions at frequencies 
which lie in each of the standard EVLA and ALMA cm to 
sub-mm observing windows. The angular scale and flux were 
corrected for a source distance of 3kpc — a typical distance 
to nearby H II regions. 

The synthetic line emission cubes at all frequencies have 
a similar morphology. Figure [5] shows the integrated inten- 
sity map of the synthetic H68a line emission cube as an 
example. At radii larger than a few arcseconds the RRL 
emission is optically thin and the profiles are symmetric. 
At smaller radii (approaching ~1") the emission intensity 
increases but starts to become optically thick so the pro- 
files begin to deviate from symmetric. At a radius of ~1" 
the emission becomes optically thick. The intensity then 
drops sharply and the profile asymmetries become more pro- 
nounced with decreasing radius. The transition from opti- 
cally thin to thick emission is seen as a very pronounced 
boundary in the RRL integrated intensity maps. 

We then simulated observing the synthetic maps with 
the EVLA and ALMA. For each telescope we chose a set 
of observational parameters (integration time, array config- 
uration, correlator setting etc.) that represent the expected 
range in these parameters for typical observations. We chose 
integration times of 0.5 hours and 8 hours to represent snap- 
shot imaging and deep integrations, respectively. We se- 
lected array configurations which gave synthesized beams 
close to either 1", 0.3" or 0.1" (corresponding to physical 
scales of 15, 4.4 and 1.5 mpc, respectively), including both 
compact and extended array configurations to investigate 



© 0000 RAS, MNRAS 000, 000-000 



Understanding hydrogen recombination line observations with ALMA and EVLA 7 



1.0 



0.9 



^ 0.8 



0.7 







-1 - - - 










1- 










10 8 








— 

- 
























ir/~~^~- 






I 


— 
- 
























































10 5 






III 


























10 4 










1 




2 


3 




4 


5 


6 


7 




8 




9 


1(L 



200 400 600 800 

v (GHz) 

Figure 5. Departure coefficients as function of frequency for transitions with An = 1 and electron number densities between 10 4 and 
10* cm -3 . The data is the same as in Figure [l] The ten ALMA bands are marked as vertical bands (as shown in Figure [4] there is no 
Q-transition in band 10 for the selected quantum numbers). The plot shows that departure from LTE conditions becomes increasingly 
relevant at the higher ALMA frequencies. 
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Figure 6. Integrated intensity map of the synthetic H68a emis- 
sion line cube (see Section |5.1| for details). The horizontal white 
line shows the cut used to create the position-velocity diagrams 
in Figure |10| The points A to E are positions at which spectra 
are extracted for Figures [7] and [9] Point E indicates the centre of 
the H II region. The colour scale shows the integrated intensity of 
the line emission initially increases towards the centre. However, 
between points B and C the emission becomes optically thick and 
the integrated intensity drops with decreasing radius. 



the effect of spatial filtering and surface brightness sensitiv- 
ity. 

For each observational setup we then used the CASA 
© 0000 RAS, MNRAS 000, 000-000 



task SIMDATA to generate synthetic visibilities (measure- 
ment sets) of the sky model. For the ALMA observations, 
thermal noise was added using the "tsys-atm" parameter, 
which models the typical atmospheric conditions above the 
ALMA site for a given precipitable water vapour (PWV) 
content. The PWV for each observing frequency was selected 
to be represenative of the conditions expected for observa- 
tions carried out at that frequency (see Table |2|. At present 
the "tsys-atm" mode only models the atmosphere above the 
ALMA site. For the EVLA observations we therefore used 
the "tsys-manual" mode taking the typical zenith opacities 
as a function of frequency from the AlP^Jtask CLCOR, as- 
suming a ground temperature of 11.7 C (the yearly median 
surface temperature at the VLA site) and a sky temperature 
of 275 We note in passing that corruption of phases us- 
ing "tsys-atm" will be worse than the errors in real ALMA 
data because the water vapour radiometers should partially 
correct for this. We also note that the above modelling does 
not take into account potential line confusion which will be 
an issue for some lines, particularly at higher frequencies. A 
sampling time of 10 seconds was used throughout for both 
EVLA and ALMA observations. 

The synthetic visibilities were inverted, cleaned and re- 
stored using the CASA task CLEAN. Dirty images were cre- 
ated with an area large enough to cover the model emission 
and with a pixel scale of ~l/3 the synthesised beam. Mode 
"MFS" with index = 1 (assuming no frequency dependence 



5 http: / / www.aips.nrao.edu / index. shtml 

6 These values were taken from the VLA test memo number 232: 
|http://www. vla.nrao.edu/memos/test/232/232.pdf 
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Table 2. Telescope, radio recombination line transitions and at- 
mospheric parameters used for the simulated observations. See 
text for details. 



H39a_ 
H68a_ 



Telescope 


Transition 


Frequency 


PWV 


'zenith 






(GHz) 


(mm) 




EVLA 


H186a 


1.013767 




0.008 


EVLA 


H117a 


4.053878 




0.01 


EVLA 


H93a 


8.045603 




0.01 


EVLA 


H68a 


20.46177 




0.05 


EVLA 


H63,3 


50.19619 




0.05 


ALMA 


H39a 


106.7374 


2.7 




ALMA 


H30a 


231.9010 


1.8 




ALMA 


H377 


346.7585 


1.3 




ALMA 


H26a 


353.6228 


1.3 




ALMA 


H24a 


447.5403 


1.3 





v lsr [km/a] 



on the flux of the continuum emission) and "channel" were 
used to clean the continuum and line cubes, respectively. 
A clean box was created based on the known structure of 
the model emission and a first-pass clean was performed on 
the cubes using a threshold of 3 times the theoretically pre- 
dicted thermal noise. These results were then inspected and, 
where necessary, optimised by additional manual cleaning. 
All the images were corrected for primary beam attenuation 
as this becomes important at the higher observing frequen- 
cies where the primary beam is comparable to the angular 
extent of the model emission. 



5.2 Results 

We find the observational results depend critically on three 
parameters: (i) whether the observations have sufficient an- 
gular resolution to resolve the optically thick central region; 
(ii) whether they have sufficient surface brightness sensitiv- 
ity to detect this optically thick line emission; (iii) whether 
the intrinsic line profiles are narrow enough to spectrally 
resolve the asymmetries. 

The issue with the intrinsic line profile being too wide to 
spectrally resolve the line asymmetries dominates at lower 
frequencies because the pressure broadening is such a steep 
function of frequency (see Equation (13 1 in Section [2J. This 
makes it impossible to detect the line asymmetries in the 
H186a observations. 

As shown in Figure [7] the RRL emission is strongly 
detected in the 1" observations. However, the line profiles 
for all transitions are nearly symmetric, even though the 
lower frequency transitions are clearly asymmetric in the 
model data cubes. This is because the angular resolution 
is insufficient to spatially resolve the line emission coming 
from the optically thin shell at large radii, from the optically 
thick emission observed towards the center of the source. 
As the optically thin emission is brighter and has a larger 
solid angle, it dominates the optically thick emission when 
they are both in the same synthesised beam. It is interesting 
to note in Figure [7] that the lower frequency H68a transi- 
tion is skewed by a few km/s to lower velocities than the 
higher frequency H39q transition, which peaks at the sys- 
temic velocity (Okm/s). The explanation for this becomes 
clear from the higher resolution observations (see below for 
more details). The H68a emission is optically thick and, due 



Figure 7. H39a and H68a spectra extracted from the central 
position ("E" in Figure [6| of the 1" angular resolution simulated 
observations (see Section j5.1| for details) . The line profiles for both 
transitions are nearly symmetric, even though the H68« line pro- 
file is clearly asymmetric in the model data cubes. This is because 
the 1" synthesised beam is not sufficient to resolve the optically 
thick central region. However, the H68o line is shifted to lower ve- 
locities than the H39« line, which peaks at the systemic velocity 
(0 km/s). This is because the H68a emission is optically thick and, 
due to the systematic velocity structure, the red-shifted emission 
is preferentially self-absorbed, leading to a blue-ward shift in the 
peak of the emission. Searching for progressive velocity offsets 
from the systemic velocity with lower frequency (more optically 
thick) RRL transitions offers a way to infer systematic motions 
in the gas, even if it is not possible to resolve the optically thick 
region. 



to the systematic velocity structure, the red-shifted emission 
is preferentially self-absorbed, leading to a blue-ward shift 
in the peak of the emission. Although the deviation from 
symmetric line profiles is only slight, searching for progres- 
sive velocity offsets from the systemic velocity with lower 
frequency (more optically thick) RRL transitions offers a 
way to infer systematic motions in the gas, even if it is not 
possible to resolve the optically thick region. However, in 
order to clearly detect the line profile asymmetries, obser- 
vations should have sufficient angular resolution to separate 
the optically thick and thin regions. Such velocity offsets 
have already been found (e.g. |Keto et al.|2008[ ). 

The 0.1" observations suffer from a different problem. 
They have sufficient angular resolution to resolve the opti- 
cally thick region but the surface brightness sensitivity in 
the 0.1" synthesised beam is too low to detect the RRL 
emission, even for an 8 hour integration with the EVLA and 
ALMA. 

The 0.3" observations have both sufficient angular res- 
olution to resolve the optically thick region and sufficient 
surface brightness sensitivity to detect the RRL emission. 
Figure|8]shows the integrated intensity (zeroth moment) and 
velocity-weighted intensity (first moment) maps of the simu- 
lated 0.3" H68q RRL observations. The integrated intensity 
map shows these observations have both sufficient surface 
brightness sensitivity to detect the RRL emission and suffi- 
cient angular resolution to resolve the optically thick region. 

Figure [9] shows the H68a, H39a and H26a spectra ex- 
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Figure 8. Integrated intensity (zeroth moment, top, with units 
Jy/beamxkm/s) and velocity-weighted intensity (first moment, 
center and bottom, with units km/s) maps of the simulated 0.3" 
H68« RRL observations. The two first moment maps have been 
made with the different intensity cuts. The central map was made 
using a threshold of 10 mjy (i.e. all pixels below 10 mjy were 
excluded before the intensity- weighted velocity was calculated). 
This accentuates the brighter emission and a clear blue-shifted 
velocity offset is seen towards the center of the region, where the 
H68a emission becomes optically thick (see Section |6] for more 
details). The bottom map was made with a much lower threshold 
of 3mjy. The lower intensity emission dominates the map. As a 
result, no velocity shift is seen. 



tracted from positions A to E in Figure [6] for the 0.3" reso- 
lution simulated observations. All spectra are plotted on the 
same scale and the systemic velocity is Okm/s. The highest 
frequency (H26«) transition is noticably brighter than the 
lower frequency transitions. This offsets the higher thermal 
noise inherent in the higher frequency observations. At the 
outer most position (A), the emission is optically thin and 
the line profiles from all transitions is symmetric and peaks 
at the same velocity. At position B the flux density from all 
transitions increases. The H39q and H26« spectra are still 
optically thin, so the line profiles are symmetric and peak 
at Okm/s. However, the H68« emission is becoming opti- 
cally thick and a red-shifted absorption shoulder is observed 
in the line profile. This effect is more pronounced at posi- 
tion C, and the H68a spectra peak velocity is blue-shifted 
to ~ —10 km/s. The flux of the higher frequency transitions 
drops compared to position B, showing they are also being 
affected by increasing optical depth. At position D the line 
profiles are all strongly self-absorbed. The measured velocity 
dispersion of the lines at this position, for example as mea- 
sured by the full width at half maximum (FWHM), would 
be much larger than at positions farther from the center. 
The self absorption in all the spectra at the central position 
(E) is strong enough to produce line profiles with two peaks 
separated by tens of km/s. Without the line profiles from the 
optically thin emission at larger radii it would not be possi- 
ble to tell if this emission was from a single, optically thick 
velocity component or two seperate, optically thin velocity 
components. 

It is important to note that the line profiles vary in a 
systematic way as the line frequency changes. At a fixed po- 
sition, the line profile becomes increasingly asymmetric as 
the frequency decreases because the optical depth increases, 
enhacing the non-LTE effects. This overcompensates the 
larger value of the departure coefficients, which are getting 
closer to the LTE- value of unity towards smaller frequencies. 
This change of the line profiles with frequency is character- 
istic for the non-LTE effect that can be used to trace sys- 
tematic gas motions. Asymmetrie^Jin the line profile due to 
asymmetries in the H II region, however, are expected to be 
equally pronounced at all frequencies. In reality, of course, 
both effects will occur simultaneously. Therefore, detailed 
radiative transfer modeling might be necessary to extract 
the desired information on systematic gas motion like infall 
or outflowj^] inside the H n region. Our implementation of 
recombination lines in RADMC-3D is an ideal tool to carry 
out such modeling. 

The information in the spectra of Figure [9] can also be 
visualised in other ways. The H68a velocity-weighted inten- 
sity (first moment) map in the central panel of Figure [8] 
shows a clear break in the velocity structure. At large radii, 
where the emission is optically thin, the intensity-weighted 

7 In principle, asymmetries could also be produced by tempera- 
ture gradients under LTE conditions. However, detailed simula- 
tions of the formation and expansion of ultracompact H II regions 
that self-consistently calculate the temperature in the ionized gas 
| |Peters et a l. 2010a c b 20lT| show that the temperature is very 
homogeneously equal to 10 4 K across the H II region. 

8 In this paper, the word outflow refers to the pressure-driven 
expansion of the H II region, not to bipolar outflows such as those 
driven from accretion processes of young stellar objects. 
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velocity is close to Okm/s — the systemic velocity of the H II 
region. Once the emission becomes optically thick, the red- 
shifted emission is preferentially absorbed due to the sys- 
tematic velocity structure (infall) of the gas. The fact there 
is more blue-shifted than red-shifted emission means the 
velocity- weighted intensity therefore shifts to ~ — lOkm/s. 
This characteristic 'bulls eye' morphology is a typical in- 
fall/outflow signature in intensity-weighted velocity maps. 
However, if the threshold used to define the integrated in- 
tensity map is dropped from lOmJy to 3mJy (a few times 
the noise level: bottom panel of Figure [9}, the intensity- 
weighted velocity is dominated by the optically thin, lower 
intensity emission at higher velocities so the 'bulls eye' pat- 
tern dissapears and a uniform velocity field at the systemic 
velocity is seen. 

The same information can be seen in the position- 



velocity (PV) diagrams in Figure 10 The H26a PV diagram 
is nearly symmetric about the systemic velocity (Okm/s). 
However, the H68a is clearly asymmetric and blue-shifted 
with respect to the systemic velocity. 
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Figure 9. H680, H39a and H26a spectra extracted from posi- 
tions A to E (top to bottom) in Figure [6] for the 0.3" resolution 
simulated observations. All spectra are plotted on the same scale 
and the systemic velocity is Okm/s. 



6 REQUIREMENTS TO OBSERVE 

ASYMMETRIC RRL PROFILES IN H 11 
REGIONS 

RRL line profile asymmetries are a potentially powerful 
probe of systematic motions (infall and outflow) in ionised 
gas. The results from the simulated observations above show 
that being able to resolve the optically thick region is cru- 
cial for any observations hoping to detect these asymmetries. 
Using a basic H n region model and making some simple as- 
sumptions we now try to predict the conditions necessary to 
observe asymmetric RRL profiles for a range of H 11 region 
properties. 

We start with a basic model of an H 11 region as a 
sphere with uniform electron number density n c , temper- 
ature T and radius R which is at a distance D. While such 
a model is undoubtedly oversimplified, these physical prop- 
erties should be thought of as similar to the beam-averaged 
values reported by observers for H 11 regions which probably 
have either multiple density and temperature components 
or temperature and density gradients. 

The critical physical scale for observing the RRL asym- 
metries is the radius at which the emission becomes opti- 
cally thick. We denote this radius i? t hick- With gas at a uni- 
form electron density and temperature, it is trivial to solve 
the one-dimensional radiative transfer equation and calcu- 
late the path length along the line of sight for which the 
emission optical depth equals unity for a given RRL fre- 
quency. We define 2 x i?thick to be the path length for which 
a photon of a given frequency, emitted from the back of the 
H 11 region would be reabsorbed before reaching the front 
edge of the H 11 region — i.e. the path length for which r = 1. 

Figure [rT] shows the dependence of i? t hick on observing 
frequency as a function of n e for H 11 regions with T = 10 4 K. 
For reference, the angular size of a source with radius i?thick 
at a distance of 2kpc is given on the right-hand vertical 
axis. Overlayed on this plot are the typical size and density 
of compact, ultracompact and hypercompact H 11 regions 
(following |Kurtz]|2005[ ), showing the regimes at which the 
emission is optically thick or thin. For the densest hyper- 

© 0000 RAS, MNRAS 000, 000-000 



Understanding hydrogen recombination line observations with ALMA and EVLA 11 





Angular Offset 



Angular Offset 



Figure 10. Position-velocity (PV) diagrams for the H26a (left) and H68a (right) cubes along the cut shown in Figure [6] for the 0.3" 
resolution simulated observations. 
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Figure 11. Threshold radius -Rthick (and corresponding angu- 
lar scale 8 at a distance of 2 kpc) at which the emission becomes 
optically thick as function of electron number denity n c and for 
frequencies of 1, 10, 100 and 1000 GHz. The crosses mark the 
regimes (from left to right) of compact, ultracompact and hyper- 
compact H II regions, respectively. 



compact H II regions, frequencies larger than 100 GHz are 
required to probe the optically thin emission. 



To distinguish the regime at which the intrinsic line 
width is too broad to spectrally resolve the line asymmetries 
we need to assume some typical velocity structure in the gas. 
The sound speed in ionized gas, c a , seems a sensible choice 
here. For example, a region undergoing pressure-driven ex- 
pansion will initially expand at the sound speed. Therefore, 
the most red-shifted and blue-shifted emission would be sep- 
arated in velocity by 2 x c a . Since the temperature inside an 
H II region is approximately T = 10 4 K everywhere, the cor- 
responding sound speed is about c s = 9km/s. As outlined in 
the Section [2j the intrinsic RRL line profile is a Voigt profile 
resulting from the contributions of a Gaussian thermally- 
broadened component and a Lorentzian pressure-broadened 
component. The width of the Voigt profile can be character- 
ized with its full width at half maximum (FWHM) / v . To 
see under which conditions the red- and blue-shifted compo- 
nents can be observed individually, /v must be compared to 
the Doppler shift due to the systematic velocity difference 
Aw, 



Av 

2 — v, 
c 



(26) 



where c is the speed of light and the factor 2 comes from 
the fact that both the red- and blue-shifted components 
have a velocity Av, but with opposite signs. This leads to 
the requirement that Av ^ /v. To be definite, we take a 
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Figure 12. Ratio of the Voigt profile FWHM fy and the Doppler 
shift Av as function of line frequency v for n c = 10 6 cm -3 and 
Ad = llkm/s. Preasure broadening dominates at low frequencies 
but becomes negligible above ~ 100 GHz. 
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Figure 13. The red- and blue-shifted line profiles and their sum 
for the same conditions as in Figure [12] at v = 100 GHz. The 
individual components can be separated. 



typical number density of a hypercompact H II region of 
n c — 10 6 cm -3 and a velocity slightly greater than the sound 
speed, An = llkm/s. Figure 12 shows the ratio fy/Av as 
function of the line frequency v. The profile width fy in- 
creases dramatically at low frequencies due to the pressure 
broadening term but approaches the thermal linewidth at 
high frequencies. The red- and blue-shifted components as 
well as their sum is shown in Figure [13] at a frequency of 
100 GHz, where the two peaks can be individually resolved. 

The figures demonstrate that pressure broadening will 
not be problematic for ALMA observations even of hyper- 
compact H n regions. On the other hand, the systematic 
gas motion must be at least as large as the sound speed 
to be detectable due to the thermal broadening. Since both 



the thermal line width (111 and the spectral separation ( 26 \ 



scale linearly with frequency v, sub-thermal motion within 
the H II region cannot be resolved independent of frequency. 
The threshold frequency where the two peaks in the line 
profile can be separated is larger the smaller the systematic 
velocity difference and the denser the H II region is. Our ex- 
ample demonstrates that even slightly super-sonic velocities 
can be separated in most ALMA bands for typical H II re- 
gion densities provided that the spectral resolution is signif- 
icantly smaller than Av. Of course, the detailed line profile 
shape and the ability for ALMA to resolve multiple veloc- 
ity components will depend on many additional factors (e.g. 
signal to noise ratio of each component, relative strength of 
each component etc.). 

Figure [14] summarises these results and those in Fig- 
ure [IT] It provides the predicted observational conditions 
required in order to detect RRL profile asymmetries from 
an H II region with T = 10 4 K. The dotted, dashed and 
solid lines show the threshold frequency for which H II re- 
gions with diameters of 0.03, 0.1 and 0.5 pc (matching the 
typical sizes of hypercompact, ultracompact and compact 
H II regions, respectively ( |Kurtz|2005[ )) at the given density 
become optically thick for free-free radiation. The angular 
scales these diameters correspond to for a source at 2kpc 
are shown in parentheses. For an H II region of a given den- 
sity and diameter, emission at frequencies below the line 
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Figure 14. Threshold frequency for the optically thick-thin tran- 
sition for H II regions with diameters 0.03, 0.1 and 0.5 pc (dotted, 
dashed and solid, respectively) and threshold frequency where 
fy = Av (dot-dash) as function of electron number density n e 
inside the H II region. Line asymmetries can be expected in the 
upper right corner of the diagram. 



will be optically thick, while emission above the line will be 
optically thin. 

The dot-dash line in Figure [14] shows the frequency at 
which fy — Av for Av = llkm/s and an H II region of a 
given density. As described above, we define this as a con- 
servative lower limit to the intrinsic linewidth at which it is 
possible to resolve line profile asymmetries. It is only at fre- 
quencies above this line, that the intrinsic linewidth would 
be small enough to detect the asymmetries. 

This important point from Figure [14] is that only the 
upper-right hand corner of the plot (above the dot-dash line 
and below the solid, dashed and dotted line) satisfies the 
conditions for the line profile asymmetries to be detected. It 
them becomes clear why strong RRL line profile asymme- 
tries have not been reported much before. The asymmetries 
have been difficult to detect due to a combination of the 
observational limits imposed by previous telescopes, in par- 
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ticular the brightness-temperature sensitivity, angular reso- 
lution and spectral resolution. 

It is clear that single-dish telescopes do not have suf- 
ficient resolution to resolve the optically thick component. 
Only interferometers pass the angular resolution criteria. At 
cm wavelengths the VLA has been the premier facility in 
terms of angular resolution and surface brightness sensitiv- 
ity. Indeed, much of the pioneering recombination line work 
was done with the VLA several decades ago at frequencies 
small er than 20 GHz (e.g. |Garay et aL|l986| |De Pree et~aT1 



1995). However, as shown in Figure 12 the pressure broaden- 



ing at these low frequencies is so large that it may have dom- 
inated the line profile and masked any asymmetries. Also, 
the limitations imposed by the previous correlator made it 
difficult to find spectral setups with simultaneously large 
enough bandwidth to fit the very broad lines (larger than 
lOOkm/s) of the high-n transitions, while at the same time 
providing the spectral resolution to resolve line asymmetries. 
Only having a few resolution elements across the line profile 
would make it difficult to identify asymmetries. These corre- 
lator restrictions have been removed with the upgrade of the 
VLA to its reincarnation as the EVLA. EVLA observations 
of RRLs towards dense H II regions at high angular/ velocity 
resolution and sensitivity may be able to detect the profile 
asymmetries. 

Higher frequency interferometers (e.g. SMA, CARMA, 
PdBI) do not suffer this band width/ velocity resolution prob- 
lem. Indeed, in going to lower n, the pressure broadening 



drops off rapidly (Equation (131, Figure 12 1. However, for 
a number of reasons the surface-brightness sensitivities of 
these interferometers is lower than that of the VLA. The 
simulated observations above show a lot of integration time 
would be required to achieve high enough resolution and, at 
the same time, sensitivity to detect asymmetric profiles. 

In summary, it is only with the order of magnitude in- 
crease in surface-brightness sensitivity, higher angular reso- 
lution and higher frequencies accessible with ALMA and the 
improved correlator with the EVLA that these may become 
routinely observable. 



7 SUMMARY 

We have presented a discussion of the capabilities for ALMA 
to observe hydrogen recombination lines and of their physi- 
cal properties. An important difference between ALMA lines 
and lines at cm-wavelengths is that the level populations 
of the former lines are farther away from LTE conditions, 
which gives rise to the formation of asymmetric line pro- 
files. We have investigated necessary conditions to observe 
such profiles with ALMA and presented synthetic ALMA 
and EVLA observations of simple model H II regions. We 
find that the asymmetric line profiles can provide useful in- 
formation about systematic gas motion, such as infall or 
outflow, within the H II region that was previously unacces- 
sible. Radiative transfer modeling, e.g. with the RADMC- 
3D implementation of recombination lines presented here, 
will be necessary to distinguish non-LTE effects from geo- 
metric effects due to asymmetries in the H II region and to 
extract the desired information from the observations. In- 
terested readers can obtain the code upon request from the 
first author. 
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APPENDIX A: 
PROFILES 



SYMMETRY OF LINE 



In this Section, we investigate under which conditions ob- 
served line profiles from symmetric H II regions are symmet- 
ric. The result will be that line profiles from isothermal, sym- 
metric H n regions under LTE conditions will always be sym- 
metric, independent of optical depth effects, whereas asym- 
metric line profiles are typical for H II regions in non-LTE 



states. This was previously noticed by Rodriguez ( 1982 1, but 



we want to make the original argument clearer. To simplify 
the calculation we will neglect the continuum opacity in this 
problem. The continuum opacity varies only mildly across 
the line profile compared to the line opacity, and we are 
interested here in true line profile asymmetries, not just a 
profile on top of an asymmetric continuum background. 

We consider a ray passing through an isothermal, spher- 
ically symmetric H II region (see Figure All. Because of 
the symmetry, all physical quantities along the ray will be 
mirror-symmetric with respect to the point where a radial 
line from the center of the H II region C cuts the ray at 
right angle. Let this point P be the origin of our coordinate 
system, and let R be the distance from P to the edge of the 
H II region along the ray, then we can write the total optical 
depth as 



a v + v> ,L 



(r) dr' . 



(Al) 




Figure Al. A ray traversing a spherically symmetric H II region. 
A radial line from the center C cuts the ray at right angle at a 
point P. The distance along the ray from P to the edge of the 
H II region is R. Because of the spherical symmetry, all physical 
quantities along the ray are mirror-symmetric with respect to the 
point P. 



and with the symmetry of the profile function (141 we find 

« I / +"'.L( _r ') = a^ ~^' : L( r ')- ( A2 ) 

By decomposing the integral ( |A1[ ) and using relation ( |A2[ ), 
we find 

t vo+v > = I {a vo+v '^{r ) + a^j->/',L(r ))dr . (A3) 
Jo 



This expression is evidently symmetric in v 1 , and hence we 
u -v' , the total optical depth of the ray is 



have t vo+v i — r, 



symmetric. This result is independent of any LTE assump- 
tions. 

The above argument demonstrates that the optical 
depth is only symmetric after integration over the whole 
region. For positions inside the H II region, say between P 



and the observer, the two contributions in Equation (A3 1 do 



not cancel because the term from the rear side is already in- 
cluded while the corresponding contribution from the front 
side is not, and a frequency dependence remains. 

In LTE, the source function of an isothermal H II region 



is simply S v = B V (T) everywhere, and thus the integral ( 25 1 
reduces to 



Iu(t v ) =fl*(T)(l-e- 



(A4) 



Since the total optical depth t v is symmetric, so is the in- 
tensity I v at the observer, neglecting the small variation of 
B V (T) over the line profile. 



In non-LTE, however, the source function (24 1 is 
S V = B V (T) h::: 



b n fin , ' 



(A5) 



where we have again neglected the continuum opacity. Thus, 
the intensity is now given by the integral 



I„(t„) = e~ T "B v {T) 



bn{Tl)Pn, m {Tl) 



dr' u . (A6) 



Because of the mirror symmetry we have v{— r') = — v(r'), 



There is no reason why this expression should still be sym- 
metric. In fact, as noted above, r v is only symmetric after 
integrating over the whole ray. The expression e r " in the 
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integral will not necessarily be symmetric, leading to a non- 
symmetric value of the integral, and hence the non-LTE de- 
parture coefficients break the symmetry of the line profil^] 
Our implementation tests (see Section B3 I demonstrate that 
this effect can be surprisingly big, proving the necessity of 
using non-LTE conditions in simulating recombination line 
observations. 



APPENDIX B: TESTS 

We verify our implementation by a series of tests conducted 



by Viner et al. (19791. All radiative transfer tests are run 



with free-free radiation and recombination lines combined. 
After the radiative transfer, we use a polynomial first-order 
fit to subtract the continuum contribution from the images. 

The tests are carried out for three different density pro- 
files, which we denote as model 1 to 3. For all three models 
we assume a constant gas temperature of T = 10 4 K and a 
microturbulent velocity of £ = 15 kms -1 . The models differ 
in the radial density profiles, which are defined as follows: 



• model 1: 

28 x 



«e 







cm" 3 0.02 pc sC r sC 0.1 pc 
otherwise, 



• model 2: 

flO 5 x exp(-700(r/pc) 2 )cm~ 3 r 0.1 pc 



1° 

model 3: 



otherwise, 



5xl0 4 cm~ 3 rsC0.24pc 







otherwise. 



(Bl) 



(B2) 



(B3) 



We also consider models with expanding and contracting 
gas velocities. These models are identical to model 1, but 
additionally include a radial gas velocity v r according to 
the following prescription: 



• model 4: 



• model 5: 



• model 6: 



v r = 80kms 1 , 
v r = —80 kms -1 , 



25 x 1 



0.1 pc 



model 7: 



-25 x 1 



0.1 pc 



kms 1 , 



kms 



(B4) 
(B5) 

(B6) 
(B7) 



There are three major differences between our implemen- 
tation in RADMC-3D and the calculations by | Viner et al.| 
(19791. First, Viner et al. (19791 used an adaptive step size 



9 Although we have neglected the continuum opacity to simplify 
the mathematics, it is clear that the effect will grow with r' v and 
thus be largest in optically thick regions. 



d 
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Figure Bl. The peak line flux relative to H85a for the recombi- 
nation lines H109a, H126a, H140a, H157a and H166a for mod- 
els 1, 2 and 3. Compare Figure 3 in | Viner et al.| l [T979| . 



to integrate the radiative transfer problem (221. In our calcu- 



lations, unless otherwise noted, we always map the model se- 
tups to a homogeneous, Cartesian grid with 100 3 grid points 



and linear dimension 0.5 pc. Second, Viner et al. ( 1979 1 de- 



termined the line intensity by comparing radiative transfer 
calculations of the full emission and absorption coefficients, 
ji/,c + jv,h and a V] c + ov,l, respectively, with reduced cal- 
culations taking only continuum contributions by j„,c and 
a„,c into account. Since this procedure is unavailable to ob- 
servers, we prefer to work with the full radiative transfer 
problem exclusively and determine the line intensities by 
continuum subtraction. Third, |Viner et al.| ( |1979[ ) include 
not only hydrogen but also helium in their models. This 
leads to a marginally larger electron number density com- 
pared to a the pure hydrogen case as in our models, and 
hence can produce slightly enhanced line intensities. 



Bl Line Strength as Function of Frequency 

To study the variation of the line strength as function of the 
transition frequency, we show in Figure[Bl]the peak line flux 
relative to the H85a line for the recombination lines H109a, 
H126q, H140Q, H157q and H166a for the models 1, 2 and 
3. These transitions span almost one order of magnitude in 
frequency, and the resulting peak line flux varies by more 
than four orders of magnitude. The line strength increases 
monotonically with decreasing principal quantum number 
number and is always the largest for model 1 and the small- 
est for model 3 for a certain transition, with model 2 lying 
in between. The small differences between Figure |B1| and 



Figure 3 in Viner et al. ( 1979 1 are likely the result of differ- 



ent methods for subtracting the continuum since the peaks 
of the higher transition lines are very weak. 



B2 Line Strength and Width as Function of Order 



Figure |B2| shows the line profiles along a ray through the 
center of the sphere of model 1 for four different recombina- 
tion lines near 5 GHz with increasing transition order An 
(H137/3, H1577, H1725 and H185e). The peak line temper- 
atures decrease and the line width increases with increasing 
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Figure B2. Line profiles for rays through the center of the sphere 
in model 1. The line temperatures arc calculated for H137/3, 
H1577, H172<5 and H185e lines. Compare Figure 2c in |Viner et al.| 
11979) . 



order of the transition. There are small deviations < 15% 



from Figure 2c in Viner et al. ( 1979 1 in the height of the 
line peaks. This could in principle be a result of insuffi- 
cient spatial resolution near the density peak in model 1 
on our homogeneous grid. To test this, we have run a radia- 
tive transfer calculation with adaptive mesh refinement for 



these lines (see Section B5 for further discussion). To resolve 
the density peak in the r _2 -profile of model 1, we introduce 
a Jeans-like refinement criterion. We require that the cell 
size da; and the electron number density n e always satisfy 
the relation 



da; 

6 x — < ( 5 x 
pc 



-1/2 



(B8) 



If this condition is not met, we split the cell under consider- 
ation into eight smaller cells with half the linear dimension. 



This procedure is iterated until equation ( B8 1 is satisfied ev- 
erywhere. This leads to the creation of five refinement levels 
on top of the 100 3 base grid. The total adaptive mesh con- 
sists of 27,282,088 cells, which is more than 27 times the 
number of cells in the homogeneous case. However, the line 
peaks are amplified only very mildly, so that the discrep- 



ancy to the results of Viner et al. ( 1979 1 remains. The same 



holds for a spherically symmetric mesh with a logarithmic 
radial coordinate axis (see Section |B5| | , so that the integra- 
tion method is unlikely to be the reason. 

Figure |B3| displays source-integrated peak line fluxes 
relative to H109a for some lines near 6 cm (H137/3, H1577, 
H172(5, H185e and H196£) as function of the transition or- 
der An for models 1, 2 and 3. The peak line flux monoton- 
ically decreases with increasing An. The peak line flux is 
the largest for model 1, the second-largest for model 2 and 
the smallest for model 3 for a fixed An. The small devia- 
tions from Figure 4a in |Viner et al.| ( |1979[ ) are most proba- 
bly introduced by the continuum subtraction since again the 
higher-order lines are very weak compared to the continuum. 

Figure |B4| shows the line width for the same source- 
integrated lines as in Figure |B3| We use the full width at 
half maximum (FWHM) as a measure of the line width. To 
determine the FWHM of the lines, we fit a Voigt profile 
to the observed line profiles and extract the Doppler width 



1.000 



"a 0.100 

as 

o 

i — i 

gf 0.010 



0.001 





Model 1 = 


- \ 


- - - - Model 2 


\V\\ 
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\VV 
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1 1 . . 





An 

Figure B3. Source-integrated peak line fluxes for H137/3, H1577, 
H172<5, H185e and H196 C relative to H109a for models 1, 2 and 
3. Compare Figure 4a in |Viner et ah] ( |1979fr . 



a and the Lorentz parameter S. The FWHM of the Voigt 
profile can be computed via the FWHM of the Gaussian 
profile /g = 2\/2 In 2a and the Lorentzian profile /l = 2S 



following Olivero & Longbothum (19771 as 



(B9) 



with the parameters Ci = 1.0692, C 2 = 0.86639 and C 3 = 
1.0. 

The line width in Figure |B4| increases monotonically 
with the order of the transition. This is because the prin- 
cipal quantum number of the transition also increases, and 



following Equation ( 13 1 pressure broadening becomes more 
important for higher quantum numbers. We also show a cal- 
culation with LTE occupation numbers and a pure Gaussian 
line profile for model 1 for comparison. While the linewidth 
for models 1 and 2 agree relatively well with Figure 4b in 
|Viner et al.| ( |19T9| >, the linewidth for model 3 is larger by a 
factor of about two in our computation for high-order transi- 



tions. This is probably because Viner et al. ( 1979 1 use a dif- 
ferent fitting formula for the 5 parameter of the Lorentzian 
profile. 

To demonstrate the goodness of the Voigt profile fit to 
the measured line profiles as well as the quality of the derived 
linewidth parameters we show the measured line tempera- 
tures for model 1 for the lines H109a (An = 1) and H196£ 
(An = 6) in Figure [B5| The line temperature for the H196£ 
is multiplied by a factor of 300 to allow a direct comparison 
to the H109q line. The Voigt fits are very good, which veri- 
fies that the derived Voigt profile parameters a and 5 have 
converged. 



B3 Line Shape as Indicator of Gas Motion 

To investigate the effect of radial gas motion on the line 
shape, we consider H85a observations of models 4 to 7. As 
explained in Appendix [A] such line profiles must always be 
symmetric with respect to the rest frequency in LTE calcu- 
lations. Figure |B6] shows the source-integrated line tempera- 
ture for the case of constant radial expansion (model 4) and 
contraction (model 5). Indeed, the LTE profiles are sym- 
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o 



— Model 1 

- Model 2 
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— LTE, Gaussian 
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An 

Figure B4. Linewidth for the same lines as in Figure |B3] and a 
control calculation of model 1 with LTE occupation numbers and 
Gaussian line profiles. Compare Figure 4b in |Viner et"aT] (1979) ■ 
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Figure B5. Source-integrated line temperatures for the H109a 
line and H196C line (multiplied by 300) for model 1. The plot 
shows both the measured profiles and the Voigt fits to the profiles. 
The fits represent the data very well. 



metric and cannot distinguish between expanding and con- 
tracting motion. The non-LTE profiles, however, are highly 



asymmetric. Viner et al. ( 1979 1 explained the asymmetric 



line profile as a result of optical depth effects (see also |Es-| 
calante et al. (1989)). As explicated in Appendix |A} it is 
crucial to carry out the calculation under non-LTE condi- 
tions to see this effect. 

Figure [B7| demonstrates the same mechanism for mod- 
els 6 and 7, this time for a line of sight through the center 
of the sphere. The strength of the effect varies with spatial 
position and is strongest for the line of sight shown. It is 
interesting that such highly asymmetric line profiles can be 
produced simply through non-LTE conditions, without any 
asymmetries in the temperature, density or velocity distri- 
bution. 



B4 Maser Amplification 

In principle, the non-LTE conditions within H II regions can 
give rise to maser amplification of the recombination line 



Model 
Model 
Model 
Model 



non-LTE 
non-LTE 
LTE 
LTE 




-100 

v (km/s) 



200 



Figure B6. Source-integrated line temperature of H85a for mod- 
els 4 and 5 in LTE and non-LTE. The LTE line profiles for expand- 
ing (model 4) and contracting (model 5) gas motion are identical 
and symmetric, whereas the non-LTE line profiles are clearly dis- 
tinguishable and asymmetric. Compare Figure 5 in |Viner et al.| 
(1979) . 
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Figure B7. Line temperature of H85a for a line of sight through 
the center of the sphere for models 6 and 7 in LTE and non-LTE. 
The LTE line profiles are symmetric and identical for model 6 
and model 7, while the non-LTE profiles differ because of their 
asymmetry. Compare Figure 6 in | Viner et aL] | |1979| ). 



emission (Gordon & Sorochenko 20021. To test the ability 
of our code to model maser amplification, we conduct the 
maser test problem of Viner et al. (19791. It uses the setup 
of model 1, but with some slight modifications. The electron 
density and temperature in the interior cavity are now n c = 
3 x 10 4 cm -3 and T = 2500 K, respectively. In the shell with 
the power-law profile, the temperature is T = 2 x f0 4 K. 
To demonstrate the maser amplification, we consider the 
emission from the shell, the cavity, and a combination of 
both separately. 

In Figure [B8l we show the source-integrated line temper- 
atures of the H85a line. The combined model that includes 
the shell and cavity simultaneously has a much larger line 
temperature than the models with the shell or cavity alone. 
As the Figure demonstrates, the line temperature of the 
combined model even significantly exceeds the arithmetic 
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Figure B8. Source-integrated line temperatures of H85a for the 
shell and cavity components as well as a combined model. The 
combined model exceeds the arithmetic sum of the components, 
demonstrating maser amplification. Compare Figure 7 in |Viner| 
|et al.|fl979) . 



sum of the line temperatures from the model components 
The agreement with Figure 7 in " 
lent. 



Viner et al. ( 1979 1 is excel 



B6 Relevance of Fine-Structure Splitting 

As already mentioned in Section [2j the recombination lines 
show fine-structure splitting for frequencies above 100 GHz, 
which is the frequency range where ALMA operates. Thus, 
it is important to know to which degree the fine-structure 
splitting affects the observational appearance of these lines. 
Although the intensity distribution of the individual 



fine-structure components is highly asymmetric ( Towle et al 



|1996p , the components are so close and the individual pro- 
files are so strongly Doppler-broadened that deviations from 
a Gaussian profile are not detectable, at least for typical H H 



region temperatures around 10 K. As Towle et al.| ( |1996 l 
show, the largest deviation of the brightest fine-structure 
component or of the intensity-weighted mean to the Ryd- 
berg frequency can be at most of the order MHz for ALMA 
frequencies. However, our numerical experiments show that 
the intensity of the superposition of all fine-structure com- 
ponents is systematically lower than the intensity given by 
the Menzel ( 1968 1 formula, and this difference can be of 



the order 10%. Since the simultaneous treatment of all fine- 
structure components is vastly more expensive than using a 



single line profile, we choose to use the Menzel ( 1968 1 ap 



proximation for ALMA frequencies but to keep in mind that 
the intensities are slightly overestimated. 



B5 Adaptive-Mesh Test 

RADMC-3D has the capability of performing radiative 
transfer calculations on adaptive meshes. This is not only 
highly demanded to post-process simulation data, but 
RADMC-3D can also perform on-the-fly refinement on se- 
tups that are entirely defined by the user, such as in our test 
cases. We have already applied this method in Section |B2| 
to better resolve the r~ 2 singularity in model 1. Here, we 
present additional tests of the adaptive mesh radiative trans- 
fer in RADMC-3D. We compare synthetic maps for model 1 
calculated on a homogeneous grid, on the Jeans-like refined 
grid presented in Section |B2[ and a grid in spherical coor- 
dinates with logarithmic radial grid spacings. The spherical 
grid has a resolution in (r, !?, (^-coordinates of 100 3 . The 
cell boundaries rj of the radial collocation points are chosen 
such that 

i/N r 

(B10) 

with the inner radius ri n = 10~ 3 pc, the outer radius r out — 
0.1 pc and the number of cells in radial direction N r = 100. 
Here, the index i runs from to N r . 

Figure |B9| displays the peak line temperature of the 
H137/3 line for model 1. As expected, the line temperature is 
the highest at the boundary of the cavity, where the electron 
number density n c has its maximum. Of course, the spherical 
shape of the shell in model 1 is much better represented in 
spherical coordinates. It is also cleary visible that the gas 
with high n e is better resolved on the adaptive mesh than 
on the homogeneous grid. Note, however, that the numerical 
values of the line temperature agree very well for the three 
different grids. 
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Figure B9. Maps of the peak line temperature of the H137/3 line for model 1. Shown are images for the homogeneous grid (left), the 
adaptivcly refined grid (middle) and the spherically symmetric grid (right). 
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